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Abstract 

According to a theorem of Poincare, the solutions to differential equations are analytic 
functions of (and therefore have Taylor expansions in) the initial conditions and various 
parameters provided that the right sides of the differential equations are analytic in the 
variables, the time, and the parameters. These Taylor expansions, which provide a transfer 
map M. between initial and final conditions, may be obtained, to any desired order, by 
integration of the complete variational equations. As an example of this approach, the 
convergence of such an expansion is investigated for the Duffing equation stroboscopic 
map in the vicinity of a infinite period doubling cascade and resulting strange attractor. 

1 Introduction 

Consider any set of m first-order differential equations of the form 

Za = fa(zi,'-- ,Zm), a =!,•••, 771. (1.1) 

Here t is the independent variable and the quantities z\, ■ ■ ■ ,z m are dependent variables. 
For notational convenience let us introduce the vector z 6 W 11 with components z±, ■ ■ ■ ,z m 
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and the vector / G M m with components /i, ■ • • , f m - With this notation, the set of equa- 
tions (1.1) can be written more compactly in the form 

z = f(z,t). (1.2) 

Let z° be a vector of initial conditions specified at some initial time t = t°, 

z(t°) = z°. (1.3) 

Then, under mild conditions imposed on the functions f a that appear on the right side of 
(1.1) and thereby define the set of differential equations, there exists a unique solution 

z a {t)= la (z°;t°,t), a = l,m (1.4) 

of (1.1) with the property 

z a (t°) = 7a (z°;t°,t ) = z° a , a = l,m. (1.5) 

Let us also introduce a vector 7 G M m with components 71, • • • , 7 m . With this notation 
(1.4) and (1.5) can be written more compactly in the form 

z(t) = 7 (z°;t°,t), (1.6) 

with 

z(t°) = 1 (z ;t ,t°) = z . (1-7) 

Now assume that / is analytic (within some domain) in z and t. By this we mean that 
the functions f a are analytic (within some domain) in the components of z and the time 
t. Then, according to a theorem of Poincare, the solution quantities z a given by (1.4) (i.e., 
the components 7 a of 7) will be analytic (again within some domain) in the components 
of the initial condition vector z° and the times t° and t [1,2,3]. In vector notation, we say 
that z is analytic in z°, t°, and t. 

Poincare established this result on a case-by-case basis as needed using Cauchy 's method 
of majorants. It is now more commonly established in general using Picard iteration, and 
appears as background material in many standard texts on ordinary differential equations 
[4,5]. It is also worth noting that Poincare's result can be generalized. Suppose, for 
example, that the functions f a are only piece-wise continuous in t but are analytic in z 
for each value of t. Then the solution will be piece-wise differentiable in t. Remarkably, 
it will remain analytic in z° [6,7]. If, piece- wise, the functions f a have some order of 
differentiability in t and analyticity in z, then the solution will have, piece-wise, one order 
more differentiability in t and analyticity in z . 

Suppose that z d (t) is some given design solution to these equations, and we wish to 
study solutions in the vicinity of this solution. That is, we wish to make expansions about 
this solution. Introduce a vector of m deviation variables £ by writing 

z = z d + (. (1.8) 
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Then the equations of motion (1.2) take the form 

z d + ( = f(z d + (,t). (1.9) 

We now assume that the right side of (1.2) is analytic about z d . Then we may write the 
relation 

f(z d + (,t) = f(z d ,t)+g(z d ,t,() (1.10) 
where each component g a of g has a Taylor expansion of the form 

g a (z d ,t,0 = ^9 r a(t)G r (C). (1.11) 

r 

We also write (1.11) more compactly in the vector form 

g(z d 1 t,C) = Y,9 r (t)G r (C). (1-12) 

r 

Here the G r (£) are the various monomials in the m components of ( labeled by an index r 
using some convenient labeling scheme, and the g r a are (generally) time-dependent coeffi- 
cients which we call forcing terms. By construction, all the monomials G r (Q occurring on 
the right sides of (1.11) and (1.12) have degree one or greater. We note that the g r a {t) are 
known once z d (t) is given. 

By assumption, z d is a solution of (1.2) and therefore satisfies the relations 

z d = f(z d ,t) (1.13) 

or, in terms of components, 

Z d a = f a (z d ,t). (1.14) 

It follows that the deviation variables satisfy the equations of motion 

Ca=9a(z d ,t,0 = ^9 r a(t)G r (C) (1.15) 

r 

or, more compactly, 

C = g(z d ,t,() = Y,9 r mr(C). (1.16) 

r 

These equations are evidently generalizations of the usual first-degree (linear) variational 
equations, and will be called the complete variational equations. 

Consider the solution to the complete variational equations with initial conditions Q 
specified at some initial time f . As described earlier, we expect that under suitable condi- 
tions this solution will be an analytic function of the initial conditions Q. Also, since the 
right side of (1.16) vanishes when all the components of £ vanish [all the monomials G r 
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in (1.16) have degree one or greater], £(i) = is a solution to (1.16). It follows that the 
solution to the complete variational equations has a Taylor expansion of the form 

Ca(t) = ^ h amr(C) (1-17) 
r 

where the h r a (t) are functions to be determined, and again all the monomials G r that occur 
have degree one or greater. When the quantities h r a {t) are evaluated at some final time 
P, (1.17) provides a representation of the transfer map Ad about the design orbit in the 
Taylor form 

( f a=Ca(t f ) = Y,K(t f )Gr(C), (1-18) 

r 

or, more compactly, 

(f = ^h r (t f )G r (C) (1.19) 

r 

where each h r (t^) is an array of m entries consisting of the functions h r a {t^) with a = 1, m. 

The organization of this paper is as follows: Section 2 treats, as an example, the Duffing 
equation and describes the properties of an associated stroboscopic transfer map Ad. These 
properties are obtained by numerical integration of the equations of motion. This work is 
necessary to get the lay of the land. For comparison, Section 3 studies some of the properties 
of the truncated Taylor maps Ad n obtained by solving the variational equations numerically. 
There we will witness the remarkable fact that a truncated Taylor map approximation to 
Ad can reproduce the infinite period-doubling Feigenbaum cascade and associated strange 
attractor exhibited by the exact Ad. Truncation at order 3 already produces a map Ad% with 
the right qualitative behavior, and agreement with exact numerical results improves as n is 
increased beyond n = 3. Thus, in effect, we make a preliminary study of some convergence 
properties of truncated Taylor approximations to the Duffing equation stroboscopic map. 
A final section provides a concluding summary and discussion. 

2 Duffing Equation Example 
2.1 Introduction 

As an example application, this section studies some aspects of the Duffing equation [8]. 
The behavior of the driven Duffing oscillator, like that of generic nonlinear systems, is 
enormously complicated. Consequently, we will be able to touch only on some of the 
highlights of this fascinating problem. 

Duffing's equation describes the behavior of a periodically driven damped nonlinear 
oscillator governed by the equation of motion 

x + ax + bx + cx 3 = dcos(Qt + tp). (2.1) 
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Here tp is an arbitrary phase factor that is often set to zero. For our purposes it is more 
convenient to set 

4> = tt/2. (2.2) 

Evidently any particular choice of ip simply results in a shift of the origin in time, and this 
shift has no physical consequence since the left side of (2.1) is independent of time. 

We assume b, c > 0, which is the case of a positive hard spring restoring force. 1 We 
make these assumptions because we want the Duffing oscillator to behave like an ordinary 
harmonic oscillator when the amplitude is small, and we want the motion to be bounded 
away from infinity when the amplitude is large. Then, by a suitable choice of time and 
length scales that introduces new variables q and r, the equation of motion can be brought 
to the form 

q + 2(3q + q + q 3 = -e sin lot, (2.3) 

where now a dot denotes d/dr and we have made use of (2.2). In this form it is evident 
that there are 3 free parameters: /3, e, and co. 



2.2 Stroboscopic Map 

While the Duffing equation is nonlinear, it does have the simplifying feature that the 
driving force is periodic with period 

T = 2tt/lo. (2.4) 
Let us convert (2.3) into a pair of first-order equations by making the definition 

p = q, (2.5) 

with the result 

q = p, 

p = —2f3p — q — q 3 — esinwr. (2-6) 

Let q°,p° denote initial conditions at r = 0, and let q ,p be the final conditions resulting 
from integrating the pair (2.6) one full period to the time r = T. Let Ai denote the transfer 
map that relates q 1 ^ 1 to q°,p°. Then, using the notation z = (q,p), we may write 

z 1 = Mz°. (2.7) 

Suppose we now integrate for a second full period to find q 2 ,p 2 . Since the right side of 
(2.6) is periodic, the rules for integrating from r = T to r = 2T are the same as the rules 
for integrating from r = to r = T. Therefore we may write 

z 2 = Mz 1 = M 2 z°, (2.8) 
Other authors consider other cases, particularly the 'double well' case b < and c > 0. 
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and in general 

z n + l = Mz 11 = M n+1 z°. (2.9) 

We may regard the quantities z n as the result of viewing the motion in the light provided 
by a stroboscope that flashes at the times 2 

r n = nT. (2.10) 

Because of the periodicity of the right side of the equations of motion, the rule for sending 
z n to z n+l over the intervals between successive flashes is always the same, namely Ai. For 
these reasons M is called a stroboscopic map. Despite the explicit time dependence in the 
equations of motion, because of periodicity we have been able to describe the long-term 
motion by the repeated application of a single fixed map. 

2.3 Feigenbaum Diagram Overview 

One way to study a map and analyze its properties, in this case the Duffing stroboscopic 
map, is to find its fixed points. When these fixed points are found, one can then display 
how they appear, move, and vanish as various parameters are varied. Such a display 
is often called a Feigenbaum diagram. This subsection will present selected Feigenbaum 
diagrams, including an infinite period doubling cascade and its associated strange attractor, 
for the stroboscopic map obtained by high-accuracy numerical integration of the equations 
of motion (2.6). They will be made by observing the behavior of fixed points as the driving 
frequency uj is varied. For simplicity, the damping parameter will be held constant at the 
value j3 = 0.10. Various sample values will be used for the driving strength e. 3 

2.3.1 A Simple Feigenbaum Diagram 

Let us begin with the case of small driving strength. When the driving strength is small, 
we know from energy considerations that the steady-state response will be small, and 
therefore the behavior of the steady-state solution will be much like that of the driven 
damped linear harmonic oscillator. That is, for small amplitude motion, the q 3 term 
in (2.3) will be negligible compared to the other terms. We also know that, because of 
damping, there will be only one steady-state solution, and therefore Ai has only one fixed 
point z* such that 

Mz f = z f . (2.11) 

(Here we use the super or subscript / to denote fixed whereas in earlier sections it denoted 
final.) Finally, again because of damping, we know that this fixed point is stable. That 

2 Note that, with the choice (2.2) for ip, the driving term described by the right side of (2.3) vanishes at 
the stroboscopic times r™. 

3 Of course, one can also make Feigenbaum diagrams in which some other parameter, say e, is varied 
while the others, including to, are held fixed. 
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is, if Ai is applied repeatedly to a point near z* , the result is a sequence of points that 
approach ever closer to z* . For this reason a stable fixed point is also called an attractor. 

Figure 1 shows the values of <z/(w) for the case e = 0.150, and Figure 2 shows Pf(oj). In 
the figures the phase-space axes are labeled as and poo to indicate that what are being 
displayed are steady-state values reached after a large number of applications of M.. As 
anticipated, we observe from Figures 1 and 2 that the response is much like the resonance 
response of the driven damped linear harmonic oscillator. 1 Note that the coefficient of q 
in (2.3) is 1, and therefore at small amplitudes, where g 3 can be neglected, the Duffing 
oscillator described by (2.3) has a natural frequency near 1. Correspondingly, Figure 1 
displays a large response when the driving frequency has the value oj — 1. Observe, 
however, that the response, while similar, is not exactly like that of the driven damped 
linear harmonic oscillator. For example, the resonance peak at oj ~ 1 is slightly tipped to 
the right, and there is also a small peak for oj ~ 1/3. 

Our strategy for further exploration will be to increase the value of the driving strength 
e, all the while observing the stroboscopic Duffing map Feigenbaum diagram as a function 
of oj. We hasten to add the disclaimer that the driven Duffing oscillator displays an 
enormously rich behavior that varies widely with the parameter values /3, e, oj, and we 
shall be able to give a brief summary of only some of it. Also, for brevity, we shall 
generally only display qf(oj). 

2.3.2 Saddle-Node (Blue-Sky) Bifurcations 

Figure 3 shows the q Feigenbaum diagram for the case of somewhat larger driving strength, 
e = 1.50. For this driving strength the resonance peak, which previously occurred at oj ~ 1, 
has shifted to a higher frequency and taken on a more complicated structure. There are 
now also noticeable resonances at lower frequencies, with the most prominent one being at 
oj ~ 1/2. 

Examination of Figure 3 shows that for oj < 1.5 there is a single stable fixed point 
whose trail is shown in black. Then, as oj is increased, a pair of fixed points is born at 
oj ~ 1.8. 5 One of them is stable. The other, whose trail as oj is varied is shown in red, is 
unstable. That is, if Ai is applied repeatedly to a point near this fixed point, the result is 
a sequence of points that move ever farther away from the fixed point. For this reason an 
unstable fixed point is also called a repellor. 

This appearance of two fixed points out of nowhere is called a saddle-node bifurcation 
or a blue-sky bifurcation, and the associated Feigenbaum diagram is then sometimes called 
a bifurcation diagram. 6 The original stable fixed point persists as oj is further increased so 

4 It was the desire for q x to exhibit a resonance-like peak as a function of ui that dictated the choice 
(2.2) for tp. 

5 Actually, in the analytic spirit of Poincare, these fixed points also exist for smaller values of u, but are 
then complex. They first become purely real, and therefore physically apparent, when w ~ 1.8. 

6 Strictly speaking, a Feigenbaum diagram displays only the trails of stable fixed points while a bifurcation 
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Figure 2: Feigenbaum diagram showing limiting values function of uj (when f3 = 0.1 

and e = .15) for the stroboscopic Duffing map. 

that over some u range there are 3 fixed points. Then, as uj is further increased, the original 
fixed point and the unstable fixed point move until they meet and annihilate when uj ~ 2. 6. 7 
This disappearance is called an inverse saddle-node or inverse blue-sky bifurcation. Finally, 
for still larger uj values there is again only one fixed point, and it is stable. 

We remark that the appearance and disappearance of stable-unstable fixed-point pairs, 
as uj is varied, has a striking dynamical consequence. Suppose, for example in the case of 
Figure 3, that the driving frequency uj is below the value uj ~ 1.8 where the saddle-node 
bifurcation occurs. Then there is only one fixed point, and it is attracting. Now suppose 
uj is slowly increased. Then, since the fixed-point solution is attracting, the solution for 
the slowly increasing uj case will remain near this solution. See the upper black trail in 
Figure 3. This "tracking" will continue until uj reaches the value uj ~ 2.6 where the inverse 
saddle-node bifurcation occurs. At this value the fixed point being followed disappears. 
Consequently, since the one remaining fixed point is also an attractor, the solution evolves 
very quickly to that of the remaining fixed point. It happens that the oscillation amplitude 
associated with this fixed point is much smaller, and therefore there appears to be a sudden 
jump in oscillation amplitude to a smaller value. Now suppose uj is slowly decreased from 
a value above the value uj ~ 2.6 where the inverse saddle-node bifurcation occurs. Then 
the solution will remain near that of the fixed point lying on the bottom black trail in 

diagram displays the trails of all fixed points. 

7 Actually, they are not destroyed, but instead become complex and therefore disappear from view. 
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Figure 3: Feigenbaum/bifurcation diagram showing limiting values qoo as a function of u 
(when f3 = 0.1 and e = 1.5) for the stroboscopic Duffing map. Trails of the stable fixed 
points are shown in black. Also shown, in red, is the trail of the unstable fixed point. 
Finally, jumps in the steady-state amplitude are illustrated by vertical dashed lines at 
U) ~ 1.8 and oj ~ 2.6. 
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Figure 3. This tracking will continue until oj reaches the value oj ~ 1.8 where the fixed 
point being followed disappears. Again, since the remaining fixed point is attracting, the 
solution will now evolve to that of the remaining fixed point. The result is a jump to 
a larger oscillation amplitude. Evidently the steady-state oscillation amplitude exhibits 
hysteresis as oj is slowly varied back and forth over an interval that begins below the value 
where the first saddle-node bifurcation occurs and ends at a value above that where the 
inverse saddle-node bifurcation occurs. 

2.3.3 Pitchfork Bifurcations 

Let us continue to increase e. Figure 4 shows that a qualitatively new feature appears 
when e is near 2.2: a bubble is formed between the major resonant peak (the one that 
has saddle- node bifurcated) and the subresonant peak immediately to its left. To explore 
the nature of this bubble, let us make e still larger, which, we anticipate, will result in 
the bubble becoming larger. Figure 5 shows the Feigenbaum diagram in the case e = 5.5. 
Now the major resonant peak and the subresonant peak have moved to larger oj values. 
Correspondingly, the bubble between them has also moved to larger oj values. Moreover, 
it is larger, yet another smaller bubble has formed, and the subresonant peak between 
them has also undergone a saddle-node bifurcation. For future use, we will call the major 
resonant peak the first or leading saddle-node bifurcation, and we will call the subresonant 
peak between the two bubbles the second saddle-node bifurcation, etc. Also, we will call 
the bubble just to the left of the first saddle- node bifurcation the first or leading bubble, 
and the next bubble will be called the second bubble, etc. 

We also note that three short trails have appeared in Figure 5 just to the right of 
oj = 4. They correspond to a period-three bifurcation followed shortly thereafter by an 
inverse bifurcation. Actually, much closer examination shows that there are six trails 
consisting of three closely-spaced pairs. Each pair comprises a stable and an unstable fixed 
point of the map M 3 . They are not fixed points of M itself, but rather are sent into each 
other in cyclic fashion under the action of M . 

Figure 6 shows the larger (leading) bubble in Figure 5 in more detail and with the 
addition of red lines indicating the trails of unstable fixed points. It reveals that the bubble 
describes the simultaneous bifurcation of a single fixed point into three fixed points. Two 
of these fixed points are stable and the third, whose q coordinate as a function of oj are 
shown as a red line, is unstable. What happens is that, as oj is increased, a single stable 
fixed point becomes a triplet of fixed points, two of which are stable and one of which is 
unstable. This is called a pitchfork bifurcation. Then, as oj is further increased, these three 
fixed points again merge, in an inverse pitchfork bifurcation, to form what is again a single 
stable fixed point. 
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Figure 5: Feigenbaum diagram showing limiting values as a function of oo (when f3 = 0.1 
and e = 5.5) for the stroboscopic Duffing map. The first bubble has grown, a second smaller 
bubble has formed to its left, and the subresonant peak between them has saddle-node 
bifurcated to become the second saddle-node bifurcation. 
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Figure 6: An enlargement of Figure 5 with the addition of red lines indicating the trails of 
unstable fixed points. 
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2.3.4 A Plethora of Bifurcations and Period Doubling Cascades 

We end our numerical study of the Duffing equation by increasing e from its earlier value 
e = 5.5 to much larger values. First we will set e = 22.125. Based on our experience so far, 
we might anticipate that the Feigenbaum diagram would become much more complicated. 
That is indeed the case. Figure 7 displays c/oo when /3 = 0.1 and e = 22.125, as a function 
of oj, for the range oj 6 (0, 12). Evidently the behavior of the attractors for the stroboscopic 
Duffing map, which is what is shown in Figure 7, is extremely complicated. There are now 
a great many fixed points both of Ai itself and various powers of Ai. For small values of oj, 
and as in Figures 3 through 5, there are many resonant peaks and numerous saddle-node 
and pitchfork bifurcations. For larger values of oj there are more complicated bifurcations. 
In this figure, and some subsequent figures, the coloring scheme is chosen to guide the eye 
in following bifurcation trees with colors changing when the period changes. Points with 
period one are colored red and points of very high or no discernible period are colored 
black. 

Of particular interest to us are the two areas around oj = .8 and oj = 1.25. They 
contain what has become of the first two bubbles in Figure 6, and are shown in greater 
magnification in Figure 8. What has happened is that bubbles have formed within bubbles, 
and bubbles have formed within these bubbles, etc. to produce a cascade. However, these 
interior bubbles are not the result of pitchfork bifurcations, but rather the result of period- 
doubling bifurcations. For example, the bifurcation that creates the first bubble at oj ~ 1.2 
is a pitchfork bifurcation. But the successive bifurcations within the bubble are period- 
doubling bifurcations. In a period-doubling bifurcation a fixed point that is initially stable 
becomes unstable as oj is increased. When this happens, simultaneously two stable fixed 
points of Ai 2 are born. They are not fixed points of At itself, but rather are sent into 
each other under the action of Ai. Hence the name "period doubling". The map Ai must 
be applied twice to send such a fixed point back into itself. In the next period doubling, 
fixed points of Ai 4 are born, etc. However we note that, as oj increases, the sequence of 
period-doubling bifurcations only occurs a finite number of times and then undoes itself. 

Remarkably, when e is just somewhat larger, infinite sequences of period doubling 
cascades can occur. Figure 9 shows what happens when e = 25. 
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Figure 7: Feigenbaum diagram showing limiting values as a function of u (when f3 
and e = 22.125) for the stroboscopic Duffing map. 
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Figure 8: Enlarged portion of the Feigenbaum diagram of Figure 7 displaying limiting 
values qoo as a function of us (when f3 = 0.1 and e = 22.125) for the stroboscopic Duffing 
map. It shows part of the first bubble at the far right, the second bubble, and part of a 
third bubble at the far left. Examine the first and second bubbles. Each initially consists 
of two stable period-one fixed points. Each also contains, as u) is increased, the beginnings 
of period-doubling cascades. These cascades do not complete, but rather cease and then 
undo themselves by successive mergings as u> is further increased. The final result is again 
a pair of stable period-one fixed points. There are also many higher-period fixed points 
and their associated cascades. 
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Figure 9: Feigenbaum diagram showing limiting values as a function of u (when f3 
and e = 25) for the stroboscopic Duffing map. 
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2.3.5 More Detailed View of Infinite Period Doubling Cascades 

To display the infinite period doubling cascades in more detail, Figure 10 shows an en- 
largement of part of Figure 9. And Figures 11 and 12 show successive enlargements of 
parts of the first bubble in Figure 9. From Figure 10 we see that the first bubble forms as 
a result of a pitchfork bifurcation just to the right of w = 1.2, and from Figures 11 and 12 
we see that the first period doubling bifurcation occurs in the vicinity of uj = 1.268. From 
Figure 12 it is evident that successive period doublings occur an infinite number of times 
to ultimately produce a chaotic region when uj exceeds uj ~ 1.29. 




Figure 10: Enlargement of a portion Figure 9 showing the first, second, and third bubbles. 
The period-doubling cascades in each of the first and second bubbles complete. Then they 
undo themselves as uj is further increased. There is no period doubling in the third bubble 
when e = 25. 

2.3.6 Strange Attractor 

As evidence that the behavior in this region is chaotic, Figures 13 and 14 show portions 
of the full phase space, the q,p plane, when uj = 1.2902. Note the evidence for fractal 
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Figure 11: Detail of part of the first bubble in Figure 10 showing upper and lower infi- 
nite period-doubling cascades. Part of the trail of the stable fixed point associated with 
the second saddle-node bifurcation accidentally appears to overlay the upper period dou- 
bling bifurcation. Finally, there are numerous cascades followed by successive mergings 
associated with higher-period fixed points. 
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Figure 12: Detail of part of the upper cascade in Figure 11 showing an infinite period- 
doubling cascade, followed by chaos, for what was initially a stable period-one fixed point. 
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structure. The points appear to lie on a strange attractor. 
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Figure 13: Limiting values of ^oo, Poo for the stroboscopic Duffing map when oo = 1.2902 
(and (3 = .1 and e = 25). They appear to lie on a strange attractor. 
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Figure 14: Enlargement of boxed portion of Figure 13 illustrating the beginning of self- 
similar fractal structure. 
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3 Polynomial Approximations to Duffing Stroboscopic Map 



In this section we will compare the behavior of the exact map M with polynomial (trun- 
cated Taylor series) maps M n obtained by integrating the complete variational equations. 
In particular, we will make this comparison in the vicinity of the infinite period doubling 
cascade displayed in Figures 12 through 14 since duplicating the intricate behavior in this 
vicinity would seem to be particularly challenging. As an additional challenge, we will 
employ polynomial maps that involve expansions in the parameter oj as well. This is done 
by introducing an additional deviation variable £3 by writing 

u = uj d + ^ (3.1) 

and augmenting the complete variational equations by the equation 

Cs = 0. (3.2) 

Here oj d is the value of oj for the design solution. By an application of Poincare's theorem 
in the case of parameter dependence, the solution to 

z = f(z,t,\), (3.3) 

where A denotes a set of parameters, is analytic in any parameter providing / is analytic in 
that parameter. The dependence of the right sides of (2.6) is manifestly analytic in ui and 
we therefore know that the solution to (2.6) will be analytic in £3. For further discussion of 
this point, and a description of how the complete variational equations can be integrated 
to obtain the desired truncated Taylor map, see the references [9,10]. 

3.1 Performance of High-Order Polynomial Approximation 

Let Ms denote the 8 th -order polynomial map (with parameter dependence) approximation 
to the stroboscopic Duffing map M. Provided the relevant phase-space region is not too 
large, we have found that Ms reproduces all the features, described in Subsection 2.3, 
of the exact map [10]. (The phase-space region must lie within the convergence domain 
of the Taylor expansion.) This reproduction might not be too surprising in the cases of 
elementary bifurcations such as saddle-node and pitchfork bifurcations. What is more 
fascinating, as we will see, is that Ms also reproduces the infinite period doubling cascade 
and associated strange attractor displayed in Figures 12 through 14. 

3.1.1 Infinite Period Doubling Cascade 

Figure 15 shows the partial Feigenbaum diagram for the map Ms in the case that /3 = 0.1 
and e = 25. The Taylor expansion is made about the point indicated by the black dot. 
This point (the initial conditions for the design solution) has the coordinates 

q hd = 1.26082, Phd = 2.05452, w bd = 1.285. (3.4) 
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It was selected to be an unstable fixed point of Ai, but that is not essential. Any nearby 
expansion point would have served as well. Note the remarkable resemblance of Figures 12 
and 15. 

We have referred to Figure 15 as a partial Feigenbaum diagram because it shows only 
qoo and not poo. In order to give a complete picture, Figure 16 displays them both. 




Figure 15: Partial Feigenbaum diagram for the map M%. The black dot marks the point 
about which M. is expanded to yield M.% 
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OJ 

Figure 16: Full Feigenbaum diagram for the map The black dot again marks the 

expansion point. 
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3.1.2 Strange Attractor 



As displayed in Figures 17 through 20 like Ai, appears to have a strange attractor. 
Note the remarkable agreement between Figures 13 and 14 for M and their Mg coun- 
terparts, Figures 17 and 18. In the case of Ai% we have been able to obtain additional 
enlargements, Figures 19 and 20, further illustrating a self-similar fractal structure. Analo- 
gous figures are more difficult to obtain for the exact map M. due to the excessive numerical 
integration time required. By contrast the map A4s, because it is a simple polynomial, is 
easy to evaluate repeatedly. 
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Figure 17: Limiting values of goo, Poo for the map when uj = 1.2902. They appear to 
lie on a strange attractor. 
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Figure 20: Enlargement of boxed portion of Figure 19 illustrating the further continuation 
of self-similar fractal structure. 
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3.2 Performance of Lower-Order Polynomial Approximations 

We close this section with illustrations of the performances of M3 and third and fifth 
order polynomial approximations (including parameter dependence) to the exact map M.. 
All expansions are made about the point (3.4). Comparison of these performances gives 
some feeling for the convergence properties of the Taylor approximation to M. 

3.2.1 Performance of M3 

Figure 21 shows the M3 counterpart to Figure 15 produced using Aig. Evidently the 
qualitative features of the period doubling cascade are the same. Also, we have found that 
there is not qualitative agreement if M2 is used. We conjecture that generically third- 
order information is necessary and sufficient to obtain qualitative agreement for a period 
doubling cascade arising from what once was a period-one fixed point. 

Note also that M.3 does not reproduce the three features near oj = 1.265 seen in Figure 
12 for the exact Ai and in Figure 15 for Ais- We have found that these features first 
appear at n = 5. They belong to what was initially a period-three fixed point for Ai. 

Figures 22 and 23 show the M.3 counterparts to Figures 17 and 18 produced using Ais- 
Evidently there is qualitative agreement. The attractors in Figures 22 and 17 look similar. 
And, when enlarged, both show evidence of fractal structure. Compare Figures 23 and 18. 

3.2.2 Performance of At*, 

Figure 24 shows the AA5 counterpart to Figure 15 produced using Ais- Now there is 
improved quantitative agreement as well as qualitative agreement. Also, there are now 
three features near 00 = 1.265 that resemble those seen in Figures 12 and 15. 

Figures 25 and 26 show the Ai§ counterparts to Figures 17 and 18 produced using Mg. 
Again there is improved quantitative agreement. 
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Figure 21: Partial Feigenbaum diagram for the map The black dot marks the point 
about which M is expanded to yield .M3 
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Figure 22: Limiting values of goo, Poo for the map M3 when uj = 1.2902. They appear to 
lie on a strange attractor. 
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Figure 23: Enlargement of boxed portion of Figure 22 illustrating the beginning of self- 
similar fractal structure. 
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Figure 24: Partial Feigenbaum diagram for the map A4§. The black dot marks the point 
about which Ai is expanded to yield .M5 
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4 Concluding Summary and Discussion 



Poincare analyticity (and its generalization) implies that transfer maps Ai arising from 
ordinary differential equations can be expanded as Taylor series in the initial conditions and 
also in whatever parameters may be present. Section 1 described the complete variational 
equations, and described how the determination of these expansions is equivalent to solving 
the complete variational equations. Section 2 provided an overview of the properties of the 
stroboscopic transfer map Ai for the Duffing equation. Section 3 described examples of how 
n th degree approximations Ai n to Ai (including parameter dependence) could reproduce 
various features of the exact Ai. In particular it illustrated, remarkably, that Ais produced 
an infinite period doubling cascade and apparent strange attractor that closely resembled 
those of the exact map. It also illustrated how the accuracy of Ai n improves with increasing 
n. 

We have seen that there are situations in which a truncated Taylor map well reproduces 
results obtained by the integration of differential equations. This is comforting since the 
behavior of polynomial maps, because such maps can easily be evaluated repeatedly, is 
often studied in detail with the hope that the behavior of such maps is illustrative of 
what can be expected for maps in general, including the maps that arise from integrating 
differential equations. 

In view of this success, one might wonder if there are situations in which the use of 
truncated Taylor maps could replace or at least complement direct numerical integration. 
There is, of course, the question of convergence for Taylor series, and the convergence 
domain is related to the (generally unknown) singularity structure of the solution to the 
differential equation in the complex domain [10]. However, if satisfactory approximation 
can be illustrated by the comparison of numerical integration results with truncated Taylor 
results for representative solutions in some domain, then the use of truncated Taylor maps 
to find additional results may be faster than continued numerical integration. 

For example, in the case of the Duffing equation, although the determination of the 
relevant h r a (t) requires the simultaneous numerical integration of a large number of differ- 
ential equations, these equations need be integrated over only one drive period. Once the 
truncated Taylor series stroboscopic map has been found, its evaluation for any phase space 
point and any parameter value is essentially free. All that is required is the evaluation of 
two n-degree polynomials (one for (( and one for (| , the deviation variables associated 
with qf and pf , respectively) in three variables (({, Q, and Q). By contrast, the direct 
construction of a Feigenbaum diagram requires the integration of the Duffing equation for a 
large number of drive periods and a large number of parameter values. And, determination 
of the strange attractor associated with the Duffing equation requires the integration of 
the Duffing equation over thousands of drive periods. 

Suppose T2 is the time required to integrate two equations over a drive period. In our 
example, it is the time required to integrate the Duffing pair of differential equations (2.6) 
over one drive period. Suppose Tjv e is the time required to integrate N e equations over 
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one drive period. Let L(m, n) be the number of monomials of degree through n in m 
variables. It is given by the binomial coefficient 



n 



,m + n\ 

L(m,n)=[ . (4.1) 



See [10]. When working with m variables through terms of degree n, the number N e of 
differential equations to be integrated to determine the relevant functions h r a (t) is given by 
the relation 

N e = mL(m,n), (4.2) 

which amounts to 

N e = 3L(3, 8) = 3 x 165 = 495 (4.3) 

in the case of M.% for the Duffing equation including parameter dependence. We have 
found in our numerical studies that there is the approximate scaling relation 

T Ne ~ (N e /2)T 2 (4.4) 

for n < 9. That is, the computation time scales with the number of equations to be 
integrated. We conclude that in this example the use of Ais becomes advantageous once 
the number of drive periods times the number of parameter values exceeds 495/2 ~ 250. 

With regard to providing complementary information, it is common practice to in- 
tegrate the first degree variational equations in order to establish the linear stability of 
solutions. Integration of the higher degree variational equations, including possible param- 
eter dependence, provides information about nonlinear behavior /stability. As examples, 
such information is required for the control of orbits in accelerators and the understanding 
and control of aberrations in optical systems. 

In conclusion, there are applications for which use of the higher degree variational 
equations is advantageous, and the whole subject of the usefulness of truncated Taylor 
maps merits continued study. 

5 Acknowledgments 

D. Kaltchev wishes to thank his colleagues from TRIUMF and CERN, especially Richard 
Abram Baartman, for their interest and support. 



References 

[1] J. Barrow-Green, Poincare and the Three Body Problem, American Mathematical 
Society (1997). 



39 



[2] F. Browder, Edit., The Mathematical Heritage of Henri Poincare, Proceedings of Sym- 
posia in Pure Mathematics of the American Mathematical Society 39, Parts 1 and 2, 
American Mathematical Society (1983). 

[3] H. Poincare, New Methods of Celestial Mechanics, Parts 1, 2, and 3. (Originally pub- 
lished as Les Methodes nouvelles de la Mechanique celeste.) American Institute of 
Physics History of Modern Physics and Astronomy, Volume 13, D. L. Goroff, Edit., 
American Institute of Physics (1993). See page 145 of the Editor's Introduction for 
a statement of Poincare's holomorphic lemma, and Sections §20 through §27 for 
Poincare's proof. 

[4] Francis J. Murray and Kenneth S. Miller, Existence Theorems for Ordinary Differential 
Equations, New York University Press and Interscience Publishing Co. (1954). 

[5] Y. Ilyashenko and S. Yakovenko, Lectures on Analytic Differential Equations, Ameri- 
can Mathematical Society (2008). See Theorem 1.1. 

[6] E. Coddington and N. Levinson, Theory of Ordinary Differential Equations, McGraw- 
Hill (1955). See Chapter 1, Section §8, Theorem 8.4. 

[7] W. Walter, Ordinary Differential Equations, Springer (1998). See Section §13. 

[8] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and 
Bifurcations of Vector Fields, Springer (1983). 

[9] D. Kaltchev and A. J. Dragt, "Poincare Analyticity and the Complete Variational 
Equations", Physica D: Nonlinear Phenomena 242 (2013). 

[10] A. J. Dragt, Lie Methods for Nonlinear Dynamics with Applications to Accelerator 
Physics, (2012), available at http://www.physics.umd.edu/dsat/ 



40 



